Quantum Classical Correspondence for a non-Hermitian Bose-Hubbard Dimer 
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We investigate the many-particle and mean-field correspondence for a non-Hermitian Af-particle Bose- 
Hubbard dimer where a complex onsite energy describes an effective decay from one of the modes. Recently 
a generalized mean-field approximation for this non-Hermitian many-particle system yielding an alternative 
complex nonlinear Schrodinger equation was introduced. Here we give details of this mean-field approximation 
and show that the resulting dynamics can be expressed in a generalized canonical form that includes a metric 
gradient flow. The interplay of nonlinearity and non-Hermiticity introduces a qualitatively new behavior to the 
mean-field dynamics: The presence of the non-Hermiticity promotes the self-trapping transition, while damp- 
ing the self-trapping oscillations, and the nonlinearity introduces a strong sensitivity to the initial conditions in 
the decay of the normalization. Here we present a complete characterization of the mean-field dynamics and 
the fixed point structure. We also investigate the full many-particle dynamics, which shows a rich variety of 
breakdown and revival as well as tunneling phenomena on top of the mean-field structure. 
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I. INTRODUCTION 

In the past decade the theoretical investigation of Bose- 
Einstein condensates led to a widespread interest in nonlinear 
quantum theories such as the nonlinear Schrodinger equation 
of Gross-Pitaevskii type 1 1]. In contrast to nonlinear general- 
izations of quantum mechanics at a fundamental level |2|], in 
the context of ultracold atoms the nonlinearity arises as the 
consequence of an effective single particle description in a 
mean-field approximation of an initially linear many-particle 
quantum system. This limit is formally similar to the classi- 
cal limit of standard single particle quantum mechanics. In 
this spirit the mean-field approximation is often formulated as 
a replacement of the particle creation and annihilation opera- 
tors with c-numbers that describe the amplitudes of the effec- 
tive single particle wave function. The time evolution is then 
governed by canonical equations of motion based on the fact 
that nonlinear as well as linear quantum dynamics can be for- 
mulated as special cases of classical canonical dynamics on 
the phase space of pure states, the projective Hilbert space. 
Thus, for Hermitian systems, the correspondence between the 
many-particle description and the mean-field approximation 
can be investigated in analogy with the usual quantum clas- 
sical correspondence for a single particle system H-d]. In 
particular the Bose-Hubbard dimer that models N bosons in 
only two modes, became a standard example many of whose 
features can be analytically understood formal- 

For both many-particle and single-particle quantum me- 
chanics, the Hamiltonian is usually demanded to be Hermi- 
tian for the description of closed systems. However, there is 
a rapidly growing interest in the use of non-Hermitian Hamil- 
tonians arising from different areas. The first is the field of 
open quantum systems where complex energies with negative 
imaginary parts are used to describe an overall probability 
decrease that models decay, transport or scattering phenom- 
ena (see, e.g., lfl7l - [22h and references therein). Although in 
most cases these non-Hermitian Hamiltonians are introduced 



heuristically, they can be derived in a mathematically satis- 
factory way starting from a system coupled to a continuum 
of states (see, e.g., (ijiH] and references cited therein). It 
is interesting to note that within the past decade a somewhat 
orthogonal motivation also generated considerable interest in 
the physics of non-Hermitian operators. This is based on the 
observation that a class of non-Hermitian Hamiltonians re- 
specting a certain antilinear symmetry, often referred to as 
-symmetry, yields purely real eigenvalues in some param- 
eter regions 112411 . Further, with the introduction of an appro- 
priate inner product they can be used to define a fully consis- 
tent quantum theory for closed systems ll25ll . The so-called 
21 -symmetric Hamiltonians have been the subject of exten- 
sive studies in the past decade see, e.g. Ii26l1 . Recently there 
is increasing interest m<P<T -symmetric systems in the context 
of optics [27-33], where first experimental results could be 
obtained |34|,|35|]. Non-Hermitian quantum dynamics differ 
drastically from their unitary counterparts, and their generic 
features are far from being fully understood. In particular, 
the investigation of the quantum classical corresrjondence for 
non-Hermitian systems is only at its beginning 



Recently, considerable attention has been paid to non- 
Hermitian extensions of the Gross-Pitaevskii equation includ- 
ing an imaginary potential, in the context of scattering and 
transport behavior of BECs ll4ll - [46ll . as well as the impli- 
cations of decay or leaking boundary conditions in partially 
open traps ll47l - [5Qll . The corresponding non-Hermitian nonlin- 
ear Schrodinger equations have been formulated in an ad hoc 
manner as a complex generalization of the mean-field descrip- 
tion in the Hermitian case. However, for a many-particle sys- 
tem the generalization of the mean-field approximation in the 
presence of a complex potential is nontrivial and intimately 
related to the semiclassical limit of non-Hermitian single par- 
ticle quantum theories. Recently, a derivation starting from 
a non-Hermitian many-particle system has been presented in 
ll5lll for an open Bose-Hubbard dimer (52, EH] described by 
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the Hamiltonian 

H = E(d\d\ — a\a2) — 2iyd\d\ + v(d\d2 -{-dial) 

+ ^(d\di -d\d2) 2 . (1) 

Here dj and at are bosonic annihilation and creation operators 
for mode j, v is the coupling constant, and c is the strength of 
the onsite interaction. For convenience we assume both v and 
c to be positive in the following. The system is opened by 
making the onsite energy of mode 1 complex. Note that the 
expectation value of the particle number N = d\d\ + d\d2 is 
conserved and the opening describes a decay of the overall 
probability encoded in the normalization of the many-particle 
wave function. A direct experimental realization of the Hamil- 
tonian (Q]) can be achieved by using ultracold bosonic atoms 
in a finite double- well trap, confined by a small tunneling bar- 
rier on one side and an approximately infinite barrier on the 
other. The parameter y can then be tuned in the experiment by 
lowering or raising the tunnel barrier. An imaginary energy 
shift Oi = Hcpq — iyN transforms this non-Hermitian Bose- 
Hubbard dimer into a model that is fp T - symmetric in the un- 
biased case (8 = 0): 

tt<p<x = (e — iy)(d\d\ — d\d2) J rv{d\d2 J rd\d\) 

+ ^(d\di - dld 2 ) 2 . (2) 

In the present article we provide a detailed description of 
the mean-field approximation for this non-Hermitian many- 
particle system introduced in lBHl . Furthermore, we show that 
the mean-field dynamics can be formulated in terms of gen- 
eralized canonical evolution equations on the classical phase 
space given by the Bloch sphere. These equations consist of a 
combination of a familiar Hamiltonian flow and an additional 
gradient flow that accounts for damping. This structure was 
recently introduced as the classical limit of non-Hermitian 
quantum theories on a flat phase space Il22ll and it is likely 
that it holds for arbitrary phase space geometries. It is closely 
related also to canonical formulations of classical dissipative 
dynamics that have been investigated in the past two decades 
Il54l - [58ll . The full many-particle dynamics can be understood 
as quantum behavior on top of the generalized classical struc- 
ture, incorporating breakdown and revival phenomena as well 
as tunneling effects. 

The article is organized as follows: In section|II]we provide 
the background of the non-Hermitian single particle two-level 
system, and introduce a renormalized Bloch representation 
for the dynamics. Further, some concepts of non-Hermitian 
quantum mechanics that are of relevance in the following are 
provided. In section (TTTJ the non-Hermitian Bose-Hubbard 
dimer is introduced as a many-boson generalization of the 
non-Hermitian two-level system. In section |lVl we review the 
generalized mean-field approximation introduced in 115 ill and 
show that it can be expressed in a canonical form of dissipa- 
tive classical mechanics suggested in 112211 . We analyze the 
resulting mean-field dynamics in detail in section |V] and com- 
pare it to the full many-particle system in section |VI] We end 
with a brief summary and an outlook. 



II. THE NON-HERMITIAN TWO-LEVEL SYSTEM 

The non-Hermitian Bose-Hubbard dimer (Q]) can be re- 
garded as an N boson generalization of a single particle two- 
level system with an imaginary energy term modeling a decay 
from one of the states, which can be described by the 2 x 2 
Hamiltonian 

H= ( £ " v 2iy _ V £ ), b,vjeRj>0. (3) 

Here the state with the lower onsite energy is assumed to be 
stable and the other one to decay with a width y. The general 
case of two decaying states differs from this model only by an 
imaginary energy offset. Despite its simplicity the system ® 
incorporates many of the generic features of non-Hermitian 
quantum mechanics and was the subject of many studies in the 
past (see, e.g., i2lll59l - [62lo . In this section we briefly review 
some features of this system and a related !PT -symmetric 
model. Furthermore, we present a less familiar representation 
of the Bloch dynamics. 

The non-Hermitian two-level system is intimately re- 
lated to a prominent W -symmetric toy-model. Applying a 
constant energy shift H — » H + iyl, that is, \\f — » \|/e^, the sys- 
tem ^ can be mapped onto the Hamiltonian 

«~-(Vv +tr ). 

which is tP *T - symmetric for 8 = 0. Introducing the discrete 
parity operator 

that interchanges the two levels and the time reversal operator 
? : i — )> — i that performs a complex conjugation, we see that 
H commutes with (P ? , whereas it commutes neither with (P 
nor with f alone. Although in the general case for 8^0 the 
Hamiltonian © is not W -symmetric, to distinguish it from 
the purely decaying system © we shall refer to it as ( P ( T- 
symmetric in the following. 

The eigenvalues of the <£ 1 -symmetric two-level system are 
given by 

X ± = ±^(8-iy) 2 + v 2 =E ± - iT±. (6) 

Thus, although the Hamiltonian is not Hermitian, for certain 
parameters it has a purely real spectrum. In fact in the un- 
biased case 8 = there is a whole region in parameter space 
|y| < |v| in which the spectrum is real. This is illustrated in 
Fig.Hl which shows the eigenvalues of H $q: as a function of 
y for 8 = and v = 1 . In the regions of purely real eigenval- 
ues all eigenvectors are simultaneous eigenvectors of the fP T - 
operator; this is often denoted as unbroken (P T -symmetry. 
The eigenvalues of the decaying system ^ are always com- 
plex with a negative imaginary part, which is degenerate for 
both eigenvalues in the regions were the fp <T -symmetric sys- 
tem has a purely real spectrum. 
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FIG. 1 : Real (left) and imaginary (right) parts of the eigenvalues © 
of the fPT -symmetric two level system in dependence on the param- 
eter y for v = 1 . 



The eigenvalues of both the 9 T -symmetric and the de- 
caying systems degenerate along lines in the parameter space 
which are specified by 



8 = and v = ±y. 



(7) 



For y=0 this reduces to the so-called diabolical point of the 
Hermitian two level system [63]. At the complex degenera- 
cies for y ^ 0, the exceptional points (EP) lEHolESlHl, the 
essence of the peculiar behavior of non-Hermitian systems be- 
comes apparent. At an EP not only the eigenvalues, but also 
the eigenvectors coincide. Thus, while the eigenvectors build 
a basis of the Hilbert space outside the EP when they coincide 
at the EP they are not sufficient to span the Hilbert space. In 
other words, along the exceptional lines d) the Hamiltonian 
is not diagonalizable but equivalent to a Jordan block. The 
occurrence of EPs can have crucial impact on the physical be- 
havior of a system (see, e.g., [1^11 [28l I^^L [65M&7h > . For the fPT- 
symmetric system (|4]) the EPs mark the border to the region 
of broken -symmetry where the eigenvalues are complex 
JH. 

In the region of unbroken fp <T -symmetry the system © 
shows a pseudo-closed behavior. This means that with the 
introduction of an appropriate inner product the time evolu- 
tion can be expressed in a unitary way. However, this should 
not be confused with the conservation of the usual probability 
as it is given by the normalization of the wave function in the 
original inner product space ||\|/|| 2 . While this is conserved 
for the time evolution in an eigenstate with real energy, this 
is in general not true for an arbitrary initial state, due to the 
nonorthogonality of the eigenf unctions. 

The dynamics of a two-level quantum system can easily 
be expressed in closed form. For a time independent Hamil- 
tonian H the Schrodinger equation i\j/ = H\\f with the initial 
condition \\f(t = 0) = \|/o is solved by \\f(t) = £/(V)Vo, where 
U(t)= exp(-iHt) is the time evolution operator. For the (P f - 
symmetric two-level system © outside the EP one finds: 



U(t) 



cos(citf) — i£ 



sin(atf) 



-IV 



sin(qy) 



-IV 



sin ( co?) 



cos(«)+i^ 



(8) 



with the complex energy £ = 8 - 



plex frequency CO = 
eigenvalue difference (0 = 

= — iv) the frequency goes to zero. In this limit the time 



iy, and accordingly the com- 
which is determined by the 
— X-). However, at the EP 




FIG. 2: (Color online) Dynamics of the fPT -symmetric © (left) 
and the decaying $3$ (right) non-Hermitian two-level system with 
8 = 0, v = 1, and different values of y (from top to bottom: y = 
0.1, 0.5, 1, 1.5) for an initial state in level 1. Shown here are the 

/i I 



absolute values of the components of the wave function | 2 (blue 



dotted line) and |\|/2 | 2 (red dashed line) as well as the total probability 
n = lYl I 2 + I ¥2 1 2 (black solid line). 



evolution operator is given by 
#ep(0 = 



1-vt 
-ivt 



-lvt 
1+v* 



(9) 



The time evolution of the normalization n = | 2 + |\j/2| 2 is 
determined by the population imbalance according to the re- 
lation 



^ = -2y(| Vl | 2 -| ¥2 | 2 ). 



(10) 



From the behavior of the (P <I -symmetric system the dynamics 
of the non-Hermitian two-level system © can be found by 
applying the time dependent transformation \\f(t) — » e -yf \|/(7). 

Figure [2] shows some examples of the dynamics for differ- 
ent non-Hermiticities y with 8 = 0, v = 1 , and for an initial 
state in level 1 . The left column shows the dynamics for the 
tP *T - symmetric system (|4]) and the right column for the de- 
caying system © for the same parameter values. It can be 
seen that for the ( P ( T -symmetric system the normalization os- 
cillates for y < v with a period that increases with increasing 
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y and diverges to infinity as y approaches the EP, y = v. The 
rate of decrease of the normalization takes its maximum value 
when the population is in the first level; growth and decrease 
rates are balanced when both levels are equally populated; and 
the growth rate is maximal when the population in the second 
level is maximal. We observe that while for small values of 
y the system performs Rabi-type oscillations between the two 
levels, the population oscillations within each level become 
parallel when the EP is approached. This nicely illustrates the 
fact that a complex term in the energy cannot be regarded as 
an overall modulation of the normalization of the system, but 
rather changes the full dynamics in a dramatic way. The os- 
cillatory behavior breaks down completely at the EP where 
the period diverges. Instead we observe an algebraic growth 
of the probability. This can be obtained analytically from 
\\f(t) = #ep(0v(0)> with m e initial state in level 1, as 



n(t) = 1 - 2vt - 



-2vV. 



(11) 



For larger values of y, the -symmetry is broken and so is 
the balance between growth and decrease - the normalization 
grows exponentially. 

For the purely decaying system ©, on the other hand, we 
observe a monotonic decrease of the normalization. The de- 
cay behavior is not exponential, which is intuitively under- 
stood by recalling that the population only decays from one 
of the levels. Therefore, the decrease is determined by the 
population of this level, which varies in time if the system is 
not in an eigenstate. It is interesting to note that in contrast to 
the ( P ( T -symmetric system, we cannot detect an obvious trace 
of the presence of the EP in the decay dynamics for the non- 
Hermitian system ©. 

The similarity of the optical wave equations in waveguide 
structures to the Schrodinger equation makes it possible to ob- 
serve the described dynamics and the $ T -related phase tran- 
sition in optical waveguide structures with gain and loss. This 
has not only been investigated theoretically |27[|29ll . but has 
recently been realized experimentally JlH l35ll. 

Although the non-Hermitian Schrodinger equation does not 
preserve the normalization, it is possible to describe the dy- 
namics of the system consistently in terms of a Bloch vector 
that stays confined to the surface of the Bloch sphere through- 
out the time evolution. For this purpose we first define the 
renormalized state vector with the components 



9; 



(12) 



For both the decaying © and the W -symmetric system © 
the dynamics are then governed by the non-Hermitian (and 
nonlinear) effective Schrodinger equation: 



d Ap! 



dt \<p2 



8-iy(l -k) 



9i 

-e + iy(l + K)y I cp 2 



, (13) 



with K = | (pi | 2 — |(p2| 2 - This dynamics by definition conserves 
the normalization |(pi | 2 + |(p2| 2 = 1- We can then define the 
components of the normalized Bloch vector in the familiar 




-0.5 

-0.5 -0.5 





FIG. 3: Effective Bloch dynamics of the non-Hermitian two level 
system (\5l with v = 1 and different values of 8 and y. The two plots 
on the top and the left plot on the bottom are for the unbiased system 
with 8 = and increasing values of y = 0.75, 1, 1.25, respectively. 
The right plot on the bottom shows the dynamics for a biased systems 
with 8 = 0.1 and y= 0.75. 



way with respect to the renormalized wave function (p: 

1 Vi¥2 + ¥i¥2 



*x =^((pl<p2 + <pi<p2) 



2\|/*\|/i+\|4\|/ 2 



s y =1(9192-9192) = - , z^-r (I 4 ) 



s z =^((pl9i-9292) 



2i \|/^\|/i+\|4\|/ 2 

1 - \|/*\|/ 2 

2\|/*\|/i+\|4\|/ 2 ' 



Using this definition we can obtain the generalized Bloch 
equations of motion from (fT3l) as 



s x = -2es y + 4ysjcJz 
s y = 2zs x — 2vs z + 4ys y s z 

s z = 2vs y -y(l-4s 2 z ). 



(15) 



Here again the normalization + Sy + s% = \ is conserved by 
construction. 

The dynamics of the renormalized quantities decouple from 
the time dependency of the normalization n = \\\f\ | 2 + |\|/2 1 2 of 
the state vector which can be obtained from the Bloch dynam- 
ics via 



-4yfe + iK for© 



-4ys z n, 



for©. 



(16) 



This allows a separate investigation of both dynamics. 

The Bloch dynamics is organized according to the fixed 
points (the stationary states), which can be obtained analyt- 
ically from the real roots of the fourth order polynomial 



167*4 + 4 ( £2 + y2 ~ T 2 )^ - £ 2 = 0, 



(17) 
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where the corresponding s y and s x values are given by s y = 
2^(1 — 4s 2 ) and the normalization condition. For every pa- 
rameter set there are only two fixed points which can be of 
different types, including sinks and sources. In general the 
type of the fixed points can be identified from the behavior 
of the surrounding vector field in a systematic manner, which 
we postpone to the discussion of the general nonlinear case in 
section IVl 

In Fig. [3] we show four examples of the Bloch dynamics, 
three for an unbiased system with 8 = and different values 
of y, and one where all parameters are nonzero. In the first plot 
(top on the left) in Fig. [3] where y < v, we observe Rabi-type 
oscillations surrounding one of the two fixed points located 
at s z = 0. However, compared to the Hermitian case the pic- 
ture is deformed. The two fixed points are not centered at 
s y = corresponding to a phase difference of zero and 71 be- 
tween the amplitudes in the two levels, respectively, but with 
increasing y they approach each other along the equator to- 
ward s y = \ and s x = 0. This is connected to the fact that the 
strict fp -symmetry which enforces both s z and s y to be zero 
in the Hermitian case, with £,y = 0, is broken for y ^ and 
first replaced by the !PT -symmetry, which only demands that 
s z — 0. At the EP y = v (shown in the right plot on top in 
Fig. the two fixed points meet and the symmetry is broken. 
For even larger values of y one of the fixed points becomes 
a sink of the dynamics, and the other a source, both located 
at s z ^ 0, that is, they belong to configurations where one of 
the levels is favored despite the symmetry of the system. This 
could be denoted as a decay-trapping. With increasing values 
of y the sink approaches the south pole of the Bloch sphere 
(corresponding to the stable level) and the source approaches 
the north pole (corresponding to the level from which the de- 
cay happens). This is due to the fact that the Bloch dynamics 
describes the mean values of the remaining part of the pop- 
ulation which moves away from the center of the decay. For 
nonvanishing 8 the system is not <P T -symmetric, and the sit- 
uation is changed. In this case the fixed points change into a 
sink and a source for arbitrary small values of y. An example 
of the non-Hermitian Bloch dynamics for 8 ^ is depicted in 
the lower right plot in Fig. [3] 



III. THE NON-HERMITIAN BOSE-HUBBARD MODEL 

The non-Hermitian Bose-Hubbard dimer (Q]) can now be 
defined as the single particle non-Hermitian two-level system 
([3]) populated with N bosons, with the bosonic particle cre- 
ation and annihilation operators a- , dj for the two levels that 
fulfill the usual bosonic commutation relations = 8jk , 

[01,02] = 0. 

A non-Hermitian many-particle Hamiltonian of the present 
type does not describe the loss of individual particles. Rather, 
it describes the decrease in time of the probability to find the 
entire many-particle ensemble in the two modes. This in- 
formation is completely encoded in the normalization of the 
many-particle wave function The expectation value of 
the particle number operator ( x I / | J /V| x I / )/( v E / | v E / ) stays constant 



in time. In other words, the "decay" is regarded as a feature of 
the state, rather than of the particles. The fact that the Hamil- 
tonian (Q]) commutes with the number operator N implies that 
the matrix representation in the Fock (particle number) ba- 
sis has a block diagonal structure, that is, it does not induce 
coupling between subspaces associated with different particle 
numbers. Therefore, in what follows we shall restrict our dis- 
cussion to these subspaces of fixed N. 

In analogy with the Bloch representation of the single par- 
ticle system, the Hamiltonian (Q]) can also be expressed in 
the form of an angular momentum system. Introducing the 
angular momentum operators L x , L y and L z according to the 
Schwinger representation 

L x — \ (d\a2 + did^) , L y = ^ (d\d2 — d\a^) , 

L z = ^(d\d\ — a 2 d2)^ (18) 

which obey the usual SU (2) commutation relation 

[lx,Ly]=iLz, (19) 

and its cyclic permutations, the Hamiltonian (Q]) can be refor- 
mulated in the form 

k = 2 (8 - iy)L z + 2vL x + 2cL 2 z - iyN . (20) 

The conservation of N appears as the conservation of L 2 — 
Y (f + l), i.e. the rotational quantum number / = N/2. 

In the standard basis of the angular momentum algebra 
|/,m), which can be defined by the relations 

L±\l,m) = v / (/^m)(/±m+ l)|/,m± 1), 
L z \l,m) = m|/,m) (21) 

with / = N/2, the Hamiltonian ik takes the form of a tridiag- 
onal (N+l) x (N+ 1) -matrix. Special features of the spec- 
trum of the present model and a corresponding & ? -symmetric 
model are discussed in [15 2l 15311 . 

In the limit of vanishing particle interaction, c = 0, the 
eigenvalue equation is solvable in closed form and the spec- 
trum consists of multiples of the single particle eigenvalues: 



K = -Wy+ (2n — N)y (e - iy) 2 + v 2 , n = 0, 1, ...,7V. 

(22) 

Thus, for 8 = at y = ±v all eigenvalues degenerate simulta- 
neously. The corresponding eigenvectors also coalesce and 
this configuration thus corresponds to a full Jordan block 
structure of the Hamiltonian, that is, an EP of higher order 
113311 , As for the single particle system the unbiased (8 = 0) 
non-Hermitian Bose-Hubbard dimer can be mapped into a 
<T - symmetric model (0 by an imaginary energy shift k = 
^pt — iy^- F° r this model the ^PT -symmetry is broken at 
the EP where all eigenvalues become complex simultaneously. 
An arbitrary small interaction strength c ^ perturbs the sys- 
tem in a manner that leads to a splitting of the EP of higher 
order into a series of EPs of second order, that is, degeneracies 
of pairs of eigenvalues and the corresponding eigenvectors. 
The interaction thereby always shrinks the region of unbroken 
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FIG. 4: Real (left) and imaginary (right) parts of the eigenvalues 
X n = E n — iY n of the Bose-Hubbard Hamiltonian as a function of 
the non-Hermiticity y for v = 1, N = 13 particles and c = 0.5 /N (top) 
and c = 0.9 /N (bottom). 

W -symmetry. In Fig. [4] we show the eigenvalues for N = 13 
particles in dependence on the non-Hermiticity y for two val- 
ues of the interaction strength. It can be seen that the region 
of purely real eigenvalues shrinks with increasing interaction 
strength. Further details concerning the spectral behavior of 
the !P <T -symmetric model © can be found in 115 311 . Some gen- 
eral aspects of (P 1 -symmetric models of Lie-algebraic type as 
the present one have been presented in 16811 . 

The many-particle dynamics can be conveniently analyzed 
in terms of the angular momentum expectation values. The 
non-Hermitian generalization of the Heisenberg equation of 
motion for an operator A (which is not explicitly time depen- 
dent) is given by GEE El 

ifi^<\|/|%) = ( V |A^-^U| ¥ ) 

= W\AMv)-i(v\M+W, (23) 
where we decomposed the Hamiltonian into Hermitian and 
anti-Hermitian parts via ik = H — if, with H = H' and f = 
P", and introduced the notation [ , ] + for the anti-commutator. 
Thus, the equation of motion for the expectation value (A) = 
<\|/|A|\|/>/<\|/|\|/> reads 

i^(A) = ([A,^])-2iAi r , (24) 

with the covariance A\ r = [A, f ] + ) — (A) (f ) . In the case of 
the Bose-Hubbard dimer ([]]), we find for the dynamics of the 
angular momentum expectation values: 

£<4> = -2t{Ly) - 2c([L y ,L z } + ) - 2y{2A? A +A?^} 

±(ly) = 2e(4)+2 C ([4,Lj + )-2v(4)-2 Y {2A^ + A^} 

±{L Z ) = 2v{Ly) - 2y{2A; , -Aj v |. (25) 

and the normalization of the many-particle wave function |¥) 
decays according to 

±m) = -2y{2{L z ) + (N)}{V\y). (26) 




FIG. 5: (Color online) The left plot shows the dynamics of the expec- 
tation values of the angular momentum operator for an initial coher- 
ent state located at the north pole of the Bloch sphere (the decaying 
level) for N = 20 particles, v = 1, y = 0.01 and g = 0.5. The right 
plot shows the corresponding decay of the survival probability (full 
black curve) and the populations of site 1 (dashed red curve) and site 
2 (dotted blue curve) 




FIG. 6: (Color online) The left plot shows the decay of the survival 
probability (full black curve) and the populations of site 1 (dashed 
red curve) and site 2 (dotted blue curve) for an initial coherent state 
located at the south pole of the Bloch sphere for N = 20 particles, 
v = 1, y= 0.01 and g = 1. The right plot shows the corresponding 
expectation value of L z . 



The many-particle angular momentum dynamics becomes 
identical to the effective B loch-equations for vanishing inter- 
action, c = 0, if the initial state is coherent, as will become 
clear later. However, to account for the particle number the 
normalization of the many-particle wave function has to be 
associated with the N-th power of the single particle wave 
function. 

To get an impression of the behavior for nonvanishing in- 
teraction strengths we show an example of the many-particle 
dynamics for a small value of y and an intermediate value of 
the interaction strength c in Fig. [5] for an initial state where all 
particles are in the decaying mode, that is, a state located at 
the north pole of the Bloch sphere. The left plot in the figure 
shows the time evolution of the angular momentum expec- 
tation value and the corresponding Bloch sphere. Similarly 
to the Hermitian case [9, 69] the Bloch vector penetrates the 
Bloch sphere throughout the time evolution. The right side of 
the figure shows the decay behavior captured by the normal- 
ization of the many-particle wave function and the population 
probability ^\hj\^)jN of the two levels. The momentary 
decay rate is proportional to the expectation value of the z- 
component of the angular momentum, that is, the population 
imbalance of the two-modes. Thus, in comparison with the 
noninteracting case (that is equivalent to the behavior of the 
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single particle system investigated in the previous section) the 
staircase behavior of the decay is slightly changed: The steps 
are not completely flat, having a negative slope for all times, 
because the L z component does not reach the stable south pole 
in the depicted time interval. This behavior becomes more 
pronounced for stronger interaction strengths, as depicted for 
an example in Fig. [6l where on the right side the L z expecta- 
tion value is shown for comparison. The breakdown behavior 
in the dynamics of the full many-particle observables can be 
understood as a many particle effect on top of the mean-field 
dynamics which stays confined to the Bloch sphere and which 
we shall introduce in the following. 

IV. THE GENERALIZED MEAN-FIELD 
APPROXIMATION AND A CANONICAL STRUCTURE 

The mean-field approximation in the Hermitian case is of- 
ten formulated in close analogy with the classical approxi- 
mation of single particle quantum mechanics. That is, op- 
erators are replaced by c-numbers and commutators by Pois- 
son brackets, and thus, the Heisenberg equations are replaced 
by Hamiltonian equations. However, this analogy was hith- 
erto of little use for non-Hermitian many-particle systems, 
as the classical limit of non-Hermitian quantum dynamics it- 
self is still far from being understood. Thus, one had to re- 
sort to alternative formulations of the mean-field approxima- 
tion. For Hermitian quantum systems the classical analo g ca n 
be derived in an elegant way using coherent states 
This method has also proven useful in the investigation of the 
quantum-classical correspondence for cold atoms in optical 
lattices described by Bose-Hubbard type Hamiltonians where 
the condensed states are equivalent to SU (M) coherent states 
[3, 6, 72]. In HI a mean-field approximation using gener- 
alized coherent states was introduced for the non-Hermitian 
Bose-Hubbard dimer (Q]). Here we provide details of this gen- 
eralized mean-field approximation and connect it to a recently 
proposed classical approximation for non-Hermitian single 
particle quantum dynamics 112211 where a generalized canoni- 
cal structure arises. Although it has only been derived for a flat 
phase space, it has been shown that the mean-field approxima- 
tion for the present model can be formulated in terms of the 
proposed generalized canonical equations of motion. From 
a practical perspective, making use of the generalized canon- 
ical structure strongly simplifies the calculation yielding the 
mean-field dynamics. This is promising for the generalization 
to larger systems involving more than two states. 

The underlying idea of the generalized mean-field approx- 
imation Dill is to describe the whole ensemble of many- 
particles by only one macroscopic wave function in the limit 
of infinite particle number. In other words, we assume that 
the particles form a condensate throughout the time evolution. 
For a two-mode system the fully condensed states can be ex- 
pressed in the form 

|jc) = (xia\ +x 2 4) iV |0,0), (27) 
with two complex coefficients that are not necessarily nor- 



malized to unity, n = \x\ | 2 + |x2| 2 - The condensed many- 
particle wave function (f2Tb is then normalized to (x\x) — 
(\x\ | 2 + |*2 1 2 )^ = n N . These states are in fact equivalent to 
the generalized SU (2) coherent states Jzll|73|], often denoted 
also as atomic coherent states. They can be constructed by an 
arbitrary SU(2) rotation R(B^) = e i8 & s ^-^ C0S< M f an ex- 
tremal Fock state ,e.g., \N), where all particles are in the first 
mode: 

\e^)=R(e^)\N). (28) 

This is equivalent to (l27l) if we set 

x\ = ^fnz _1( ^cos^, X2 = v // ^ sm §- (29) 

Thus, the mean- field approximation is equivalent to the as- 
sumption that the many-particle state, initially chosen as a co- 
herent state, remains coherent for all times of interest. This 
assumption is in fact exact if the Hamiltonian is a linear super- 
position of the generators of the dynamical symmetry group 
117 ill , in our case for vanishing interaction c = 0. This can be 
seen by calculating the action of the time evolution operator 
on an initially coherent state. For nonvanishing interaction 
it is in general an approximation yielding the mean-field dy- 
namics. The mean-field equations of motion can thus be ob- 
tained from the quantum dynamics by replacing all expecta- 
tion values with their values in coherent states and identifying 
these with the mean-field quantities. The resulting mean-field 
dynamics can be interpreted as a special case of constrained 
quantum motion B74I1 where the constraint is that the many- 
particle state is coherent. 

Let us now derive the mean-field Bloch dynamics from the 
equations of motion for the many-particle angular momen- 
tum expectation values (f25l) using the SU (2) coherent state 
approximation. The expectation values of the Li, i = x,y,z in 
terms of the coherent state coordinates x\ and X2 read: 



(x\L x \x) 


N X\X2 


(x\x) 


2 x\x\ +X\X2 


(x\L y \x) 


N X\X2 - X\x\ 


(x\x) 


2i x\x\ +X^2 ' 


(x\L z \x) 


N x\x\ — x\x2 


(x\x) 


2 X ^ X \ X^X2 



(30) 



We can identify these quantities with the components of the 
corresponding renormalized mean-field Bloch vector: 

sj = (Lj)/N. (31) 

Comparison with the definition of the mean-field Bloch vector 
in the single particle case (ITU) reveals that the coordinates of 
the coherent state can naturally be associated with the com- 
ponents of the effective single particle wave function \\f. To 
perform the mean-field approximation we further need the ex- 
pectation values of the anti-commutators appearing in (f25l) for 
SU (2) coherent states which factorize as 

([UM+) = 2(1 +5y|, 

([Li,N} + ) = 2N{U), (32) 



8 



with N = (N) . Inserting these expressions into (l25t and taking 
the macroscopic limit N —> °° with Nc — g fixed we obtain the 
desired non-Hermitian mean-field evolution equations: 



+2es x +4gs x s z 



+4ys x s z , 
-2vs z +4ys y s z , 
+2vs y -y(l-4s 2 ). 



(33) 



These nonlinear non-Hermitian Bloch equations are real val- 
ued and conserve s 2 = s 2 +s 2 +s% = 1/4, i.e. the dynamics are 
regular and confined to the Bloch sphere. The total probability 
n decays as 



n = -2y(2s z +l)n. 



(34) 



In the limit g = in which the assumption that the many- 
particle state stays coherent in time is exactly fulfilled, these 
equations reduce to the equations for the linear single parti- 
cle two level system ([T5]) . Thus, as mentioned before, this 
captures the exact many-particle dynamics in this limit. Gen- 
eralized Bloch equations related to (l33l) also appear in a dif- 
ferent context, where the influence of decoherence is inves- 
tigated 1175147911 . It should further be noted that they can be 
considered a special case of the celebrated Landau-Lifshitz 
equations with Gilbert damping appearing frequently in mag- 
netization dynamics. 

Let us now express the mean-field dynamics in the form 
of a generalized nonlinear Schrodinger equation. In terms of 
the components \\fj of the unnormalized wave function (asso- 
ciated with the coordinates Xj of the many-particle coherent 
state) this can be formulated as: 



dt y¥2 



with 



; d /VA _ /^ + gK-2iy 



lYi 



¥i 

-E-gKJ V\|/ 2 



l¥2| 2 



l¥i| 2 - 



>2 



12- 



(35) 



(36) 



The equation of motion for the normalization in this formula- 
tion is given by h = — 2y(l — K)n. While in the limit y — >• 
the wave function stays normalized and the equations are thus 
equivalent to the usual discrete nonlinear Schrodinger equa- 
tion of Gross-Pitaevskii type, the nonlinear term gets mod- 
ified due to the non-Hermiticity. Alternatively we can ex- 
press the dynamics in terms of the renormalized wave function 
q>7 = Vy/Vw: 



dt 



Ye+gK-iy(l-K) v Wcpi 

I v -8-gK+iy(l+KW Wp 2 



(37) 

with K = | (pi | 2 — |(p2| 2 - This dynamics by definition conserves 
the normalization | |(p| | 2 = 1. 

Note that the dynamics induced by the nonlinear non- 
Hermitian Schrodinger equation (l35t differs fundamentally 
from the dynamics of a discrete Gross-Pitaevskii equation 
with an additional imaginary on-site energy, where the nonlin- 
earity is determined by K = \\\f\ \ 2 — |\|/2| 2 . This latter type of 
non-Hermitian nonlinear Schrodinger equations has attracted 



considerable attention in the context of the description of scat- 
tering phenomena and the influence of leakin g bo undaries 
for Bose-Einstein condensates recently lElMTsl 55 [H. 
Furthermore, these ad hoc nonlinear non-Hermitian equations 
also appear for absorbing nonlinear waveguides W27I4291 18211 . 

In 11221 it has been shown that the mean-field approxima- 
tion of the non-Hermitian Bose-Hubbard dimer can also be 
expressed in terms of a generalized canonical structure, as we 
will review in what follows. The generalized canonical equa- 
tions of motion proposed in 112211 are of the form 



q ) = n ] vh-g 'vr 



(38) 



where p and q are canonical phase space variables and V de- 
notes the phase space gradient, £1 is the symplectic matrix 



CI - 



-1 

1 



(39) 



and G is the corresponding Kahler metric (83|, [841] on the 
relevant phase space. The classical Hamiltonian function 
9i = H — ir is given by the expectation value of the quan- 
tum Hamiltonian in the relevant coherent states. The dy- 
namics of the normalization of the original wave function 
n= \\\f\ | 2 + |\|/2| 2 is governed by the equation of motion 

n = -2Tn. (40) 

The dynamical equation (f38l) is a combination of a canonical 
symplectic flow generated by the real part H of the Hamilto- 
nian function and a canonical gradient flow generated by the 
imaginary part T. The symplectic part evidently gives rise 
to the familiar Hamiltonian dynamics of classical mechanics. 
The gradient vector with a negative sign points in the direction 
of the steepest descent of the function T. Thus, this part of the 
dynamics aims to drive the system toward the minimum of T 
and can naturally be associated with a damping. 

The generalized canonical structure can be used to directly 
calculate the mean-field dynamics without evaluating the gen- 
eralized Heisenberg equations of motion and performing the 
coherent state approximation as follows: Our classical phase 
space is given by the Bloch sphere and can be parametrized by 
the canonical variables p and q that are related to the classical 
Bloch vector via 



$x = ^V^-p 2 cos(2q) 
s y \\]\ -p 2 sin(2q) 



(41) 



We can express the expectation value of the many-particle 
Hamiltonian (OQ) in SU (2) coherent states in the variables p,q 
to find the classical Hamiltonian function: 



H = £p + vy/l - p 2 cos(2q) + ^p 2 and T = yp. (42) 

The Kahler metric on the Bloch sphere in the variables q,p is 
given by Ii22ll 



2(1 -p 






1 



(43) 
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Evaluating the generalized canonical equations 
yields 



E + gp-v- 



:COS(2<?) 



of motion 



(44) 



p = -2y(l-p 2 ) + 2vVl-p 2 sm(2q), (45) 

which is equivalent to the nonlinear Bloch equations (l33l) . 
Similar equations also appear in a related model where a dif- 
ferent mean-field approximation is applied 118511 . 

We note here that the expressions arising for the 21- 
symmetric version of the Hamiltonian Q differ from the 
present ones by a complex energy shift. Thus, since the gen- 
eralized canonical equations of motion are invariant under a 
constant energy shift (as are the usual canonical equations of 
motion of Hamilton type), the effective dynamics resulting 
from the Hamiltonian functions related to 0]) and ©, respec- 
tively, are identical, in agreement with the previous observa- 
tions. 

Note also that the nonlinear Schrodinger equation (f37l) 
can be directly formulated as generalized complex canonical 
equations of motion for the coordinates cpi , cp^, cp2, cp|. How- 
ever, here one has to take care of the constraints confining the 
dynamics to the Bloch sphere explicitly and the expression for 
the metric gets more elaborate (see Appendix lAl). Thus, for 
practical purposes the formulation in real canonical variables 
/?, q is more convenient. 



V. MEAN-FIELD DYNAMICS AND FIXED POINT 
STRUCTURE 

In this section we analyze the mean-field behavior arising 
from the interplay of non-Hermiticity and nonlinearity. The 
mean-field dynamics is organized according to fixed points, 
which correspond to stationary solutions of the nonlinear 
complex Schrodinger equation (l35l) . In contrast to the widely 
investigated behavior of vector fields in R 2 , the general fea- 
tures of vector fields on the sphere have rarely been investi- 
gated in detail. Only recently some interest in polynomial vec- 
tor fields as the present one on the two- sphere S 2 has emerged 
in the mathematical literature [86-89]. In this context it was 
shown that the upper bound of the number of fixed points for 
a general polynomial vector field of degree 2 on the sphere is 
equal to 6. 

In the present case there are at most four fixed points that 
can be obtained analytically as the roots of a fourth order poly- 
nomial similar to the Hermitian case 19011 . To see this we have 
to study the fixed point equation defined by (f33l) with s = 0, 
which provides 



(46) 



vs y = 2 1 {\-s 2 z ). 

Using this and the normalization condition s 2 + s y + s 2 = ^ 
shows that the s z coordinates of the fixed points are given by 
the real roots of the fourth order polynomial 

4(/+y 2 )4+ 4 ^ 3 +( £2 + v2 -/-T 2 K 2 -^- £2 /4=o. 

(47) 



TABLE I: Classification of fixed points according to the eigenvalues 
X\ and X 2 °f the Jacobi matrix. 







Stable node (sink) 


Index +1 




X\ , X 2 > 


Unstable node (source) 


Index +1 




X\X 2 < 


Saddle point 


Index —1 


A,i >2 = a±ip 


a<0 


Stable focus (sink) 


Index +1 




a>0 


Unstable focus (source) 


Index +1 




a = 


Center 


Index +1 



In the following we will restrict the discussion to the unbiased 
case £ = where the polynomial (l47l) becomes biquadratic 
and the fixed points are easily found analytically. The analy- 
sis can in principle be extended to the case 8 ^ in a straight- 
forward manner. For 8 = the polynomial (l47l) has the four 

solutions s 7 



0,0. 



The corresponding val- 



ues of s x and s y are then given by (l46l) and the normalization 
condition. In summary, this yields the solutions 



V 



1 

X 

2v 





V 



Sf± = 



) 



2{ g 2 +f) 



■ 2 +f) 



(48) 



Since the components of the Bloch vector are by definition 
real valued, only the real solutions correspond to actual fixed 
points. Due to the non-Hermiticity these are not necessarily 
elliptic fixed points or saddle points, which are the only pos- 
sibilities in Hamiltonian systems. Rather, as we already ob- 
served for the linear non-Hermitian case, the additional gra- 
dient flow can lead to a destruction of periodic motion and 
introduce sinks and sources to the dynamics. In principle it 
can also lead to the emergence of limit cycles |9l|,[92|] which, 
however, have not been observed in the study of the present 
system. 

For a flow on a two dimensional surface (as in the case of 
the Bloch sphere), information on the type of fixed points can 
be obtained by the surrounding linearized fields (apart from 
special cases at parameter values for which bifurcations occur, 
see e.g., i9lll and references therein). We will now briefly in- 
troduce the classification scheme; further details can be found, 
e.g., in JH-il. 

Suppose we have a system of two first order differential 
equations which can be written in the form 



q\ =Fi(q h q 2 ), qi = F 2 (q\,q 2 ). 



(49) 



The linearization of this system around an arbitrary point is 
determined by the Jacobi matrix Dij = dFi/dqj, (ij= 1,2) 
of the vector field F at that point. The eigenvalues X\ i2 of this 
matrix at a fixed point of the dynamics, that is, a singular point 
of the vector field F, can yield information about the fixed 
point type. These eigenvalues are either real or form a com- 
plex conjugate pair, due to the reality of the matrix. One can 
distinguish four basic fixed point types (nodes, saddle points, 
foci and centers) and subclasses according to the values of 
^12, which are summarized in Table H 
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FIG. 7: Parameter regions belonging to different fixed point config- 
urations of the non-Hermitian mean-field dynamics ( [33] ). 

The so-called (Poincare) index, also listed in the table, is a 
further characteristic quantity of a singular point of a vector 
field with respect to an oriented surface (see, e.g., ||93l|96|] for 
details). It is defined as the number of revolutions of the vector 
field in traversing an arbitrary curve encircling the (isolated) 
singular point (and no other singular point). The index of a 
saddle point is —1, whereas the indices of nodes, foci and 
centers are all equal to + 1 . The number and type of singular 
points of a vector field and the possible bifurcation scenarios 
for a given manifold are restricted by the index theorem. It 
states that the sum of the indices of the singular points of a 
vector field on a manifold is independent of the choice of the 
vector field and equals the Euler characteristic %E, which is 
%e = 2 in the case of a sphere. 

The fixed points of our nonlinear non-Hermitian system 
(l33t can be categorized completely according to the above 
scheme. In summary, we can distinguish three regions in pa- 
rameter space, which are sketched in Fig. [71 

1. In region 1, for y 2 -\- g 2 < v 2 , we have only two fixed 
points s c ± which are located at the equator. For g ^ 
one of them is a sink, the other one is a source. They 
degenerate to centers for g = 0. 

2. In region 2, for y 2 + g 2 > v 2 and |y| < |v|, there are four 
coexisting fixed points, namely, a sink and a source, a 
center, and a saddle point. On the line y — the sink and 
source become centers, corresponding to the Hermitian 
self-trapping states. 

3. In region 3, which is defined by |y| > |v|, only the fixed 
points Sf± exist, namely a sink and a source. For posi- 
tive y we have a source on the northern hemisphere and 
a sink on the southern. In general (g ^ 0) they are foci, 
which become nodes in the linear limit. 

At the boundaries of these regions, at the critical parameter 
values, bifurcations that necessarily respect the index theorem 
occur. 




FIG. 8: Mean-field dynamics on the Bloch sphere for e = 0, v = 1 
and different values of y and g. (Left to right and top to bottom: 
y = 0.7, g = 0.7; y = 0.75, g = 3; y = 0, g = 3; y = 1.25, g = 3) 



Figure [8] shows examples of the Bloch dynamics (l33l) in 
the three different regions and on the Hermitian line. In the 
first plot (left on the top) the dynamics is shown for y = 0.7 
and g = 0.7, that is, in region 1. We observe deformed Bloch 
oscillations surrounding the two centers (index +1). If the 
non-Hermiticity is increased, the centers approach each other 
along the equator. However, before they meet one of them 
bifurcates at the critical circle g 2 + y 2 = v 2 (and y ^ 0) into a 
saddle (index —1) and two foci (index +1), s/±, one stable (a 
sink) and one unstable (a source). 

The second plot (right on the top) in Fig. [5] shows the result- 
ing dynamics above this bifurcation, however still in region 
2. Here we observe four fixed points resulting in a mixed dy- 
namics, where besides the periodic motion surrounding the re- 
maining center s c ± there are flows from the source to the sink. 
The appearance of the two fixed points Sf± can be viewed as a 
non-Hermitian self-trapping dynamics, which collapses to the 
Hermitian case in the singular limit y = 0. This is depicted 
in the third plot (left in the lower panel). Here the foci are 
replaced by centers. The Hermitian self-trapping effect arises 
as a bifurcation of the center (index +1) into a saddle (index 
— 1) and the additional centers (i ndex + 1) at the critical circle. 
Thus, the critical value g cr i t = \/v 2 — y 2 is decreased for y ^ 
compared to the Hermitian case. In other words, the presence 
of the non-Hermiticity promotes the self-trapping effect, how- 
ever, the resulting self-trapping oscillations are damped due to 
the non-Hermiticity, as we shall discuss later. 

In region 2 the dynamics is mainly organized by the sta- 
ble and unstable manifolds of the saddle point, as shown in 
Fig. [9] for g = 3 and two values of y. In the Hermitian case, 
y = 0, these manifolds form a single figure-eight curve, a sepa- 
ratrix, encircling the self-trapping regions around the two cen- 
tres Sf±. In addition, there is a third center localized at the 
equator opposite to the saddle point. 

In the non-Hermitian case the two self-trapping centers 
change into a sink close to the north-pole and a source close 
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FIG. 9: (Color online) Stable and unstable manifolds of the saddle 
point s c - for g = 3 and a Hermitian y = (left) and a non-Hermitian 
case y = 0.75 (right). The four fixed points are marked by red dots. 



to the south-pole. The saddle-point and the center at the equa- 
tor survive and the separatrix through the saddle point trans- 
forms into a single curve emanating from the source, passing 
through the saddle-point, encircling the center, passing again 
through the saddle and, finally, spiralling into the sink, as 
shown in the left panel of Fig. [9] for y = 0.75. The surface 
is divided into two regions, an area A c of oscillatory motion 
encircling the center, and the rest, the basin of attraction of the 
sink. With increasing interaction g, the area A c shrinks into a 
thin region close to the equator (note that the positions of the 
center and the saddle point are independent of g). Decreasing 
the interaction g, the sink and the source approach the saddle- 
point and meet at the critical value g cr i t . During this process, 
the area A c grows until it covers the whole sphere at g cr i t . 

For increasing y, starting from parameter region 2, the sad- 
dle point (index — 1) and the center (index +1) on the equator 
approach each other along the equator until they meet and an- 
nihilate for y = v at s = (0, 1 /2, 0) . For larger values of y, that 
is, in region 3 only the source and the sink remain, and the 
dynamics is fully governed by the flow from the former to the 
latter, as illustrated in the last plot (right in the lower panel) 
in Fig. [8] For g = 0, the transition occurs directly between re- 
gion 1 and 3 in a non-generic bifurcation at y = ±v (the EP), 
which is depicted in Fig. [3] where the two centers meet and 
simultaneously change into a sink and a source. 

In the Hermitian case, for g > g cr i t = v, we find self-trapping 
oscillations in the vicinity of the fixed points s/± . For y ^ 
these fixed points change into a sink and a source of the dy- 
namics which results in a damping of the self-trapping oscil- 
lations. Figure [TO1 illustrates this damping effect. Here we 
plot in false colors the time dependence of s z , the population 
imbalance between the two levels, as a function of the non- 
linearity g for an initial state at the south pole of the Bloch 
sphere for four different values of y. The first plot on the left 
shows the behavior in the case y = 0. We observe two dis- 
tinct regimes: For g < g sep = 2, the starting point, and hence 
the whole trajectory, is inside the area A c and the motion s z (t) 
shows a large amplitude oscillation extending to the vicinity 
of the north pole. For g > g sep = 2. in the self-trapping re- 
gion, the motion is confined to the neighborhood of the south 
pole. At g S ep = 2 the separatrix passes through the south pole 
and the motion starting there approaches in infinite time the 
saddle point along the stable manifold. (Note the increase 




FIG. 10: (Color online) Mean-field dynamics of s z plotted in false 
colors in dependence on the nonlinearity g for v = 1, 8 = where 
the initial state is the south pole of the Bloch sphere, for different 
values of the non-Hermiticity (top to bottom and left to right: y = 0, 
y= 0.2 < v, y= 0.75 < v, andy= 1.1 > v). 



of the period of oscillation for g — » g sep where the period di- 
verges.) This behavior continues for y ^ 0, however with a 
smaller value of g sep . For small nonlinearities the population 
is completely transferred between the two levels, that is, the 
Bloch vector oscillates between the south and the north pole 
and above g sep the oscillation stays closer to the south pole 
with increasing interaction. As observed in Fig. [U the self- 
trapping states are then associated with a sink and a source 
of the dynamics. Therefore, for a nonvanishing but subcrit- 
ical non-Hermiticity < y < v, the system relaxes to a state 
with excess population in the non-decaying state above a crit- 
ical value g sep of the interaction. This appears as a damping 
in the self-trapping oscillations, which is visible in the second 
and third plot (top right and bottom left) in Fig.[l0l A similar 
observation was reported in [97] where the effect was related 
to decoherence. For even larger values y > v in region 3, as 
shown in the right plot in the lower panel in Fig. [TUJ the os- 
cillatory motion is already destroyed even in the linear case 
g = 0, and the dynamics is dominated by the flow from the 
sink to the source irrespective of the nonlinearity, that is, the 
system stays confined to the lower half of the Bloch sphere. 

For the non-Hermitian system the normalization n, which 
can be interpreted as the "survival probability" of the system, 
is also time dependent. For the non-Hermitian two level sys- 
tem (l35t the dynamics are governed by the equation of motion 
(I34L that is, n = —2y(2s z + l)n, which does not explicitly de- 
pend on the nonlinearity. Yet, the instantaneous decay rate 
is determined by the s z component of the renormalized Bloch 
vector, whose dynamics are sensitively influenced by the non- 
linear term in the Schrodinger equation. This is illustrated in 
Fig. [TT] which shows the half life time as a falsecolor plot, as 
a function of the initial conditions (p,q) for a weak decay, 
y = 0.1, and different nonlinearities. It is clearly visible that 
the nonlinearity can stabilize the system significantly for cer- 
tain initial conditions. (Note the different colorscales.) 

To understand this behavior in more detail, we investigate 
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FIG. 1 1 : (Color online) Half life time in dependence on the initial 
condition for the non-Hermitian mean-field dynamics for 8 = for 
v = 1 and y= 0.1 and different values of the nonlinearity, (from left 
to right and top to bottom, g = 0, 0.1, 0.5, 1, 1.5, 2). 



the full time evolution of the normalization for some exam- 
ples. In Fig. [12] we show the normalization n(t) = |\|/i| 2 + 
|\|/2| 2 of the wave function as a function of time for a small 
non-Hermiticity y = 0.1 and a supercritical nonlinearity g = 3 
(blue lines), in comparison to the linear evolution for g = 
(black lines). The left plot shows the dynamics for an initial 
state at the north pole of the Bloch sphere, and the right plot 
corresponds to a state initially at the south pole. We observe 
that for an initial condition at the north pole (corresponding 
to the decaying level) the decrease of the normalization is 
slightly faster due to the nonlinearity, although from time t « 5 
onward it slows down considerably. In the limit t —> °° the 
decrease becomes exponential with a very small decay coef- 
ficient. The modulations present in the linear case are much 
less pronounced here from the very beginning. Despite these 
differences, the overall decay time characterized, e.g., by the 
half life time, is not drastically changed here. If we now turn 
to the right plot and compare the nonlinear decay behavior to 
the linear one for an initial condition in the south pole of the 
Bloch sphere (corresponding to the stable level) the induced 
changes become much more pronounced. In fact, the decay 
is considerably slowed down by the nonlinearity. For longer 
times the modulations nearly vanish and the decay becomes 
approximately exponential with the same decay coefficient as 
for the initial condition in the north pole. 

This behavior can be understood in terms of the Bloch dy- 
namics discussed before. For large nonlinearities the source 
of the Bloch dynamics moves toward the north pole, which is 
connected with the decaying level and thus a sink for the prob- 
ability. The sink of the Bloch dynamics, on the other hand, 
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FIG. 12: (Color online) Decay of the mean-field normalization n 
in dependence on the nonlinearity for a small decay y= 0.1. The 
evolution for a nonlinearity of g = 3 (blue solid curves) is compared 
to the linear case (black dashed curves). The left figure shows the 
time evolution of n(t) starting in the north pole, the right figure shows 
the same for an initial state at the south pole. 



moves close to the south pole which corresponds to the sta- 
ble level. Thus, if we start the system at the south pole (right 
plot in Fig. [12]), then due to the nonlinearity it stays on the 
southern hemisphere (— \ < s z < 0) and spirals into the sink 
of the dynamics instead of performing Rabi oscillations ex- 
tending over all values of s z . Hence the instantaneous decay 
rate is smaller than for the linear case and the decay is signif- 
icantly decelerated. If we start the dynamics at the north pole 
(left plot in Fig. [T2l) . on the other hand, the Bloch dynamics 
also move toward the sink close to the south pole, where they 
remain. However, until the small instantaneous values of the 
decay coefficient associated with the southern hemisphere of 
the Bloch sphere become relevant, the normalization already 
decayed considerably. 

Summarizing, the interplay of nonlinearity and non- 
Hermiticity introduces a qualitatively new behavior to the 
mean-field dynamics. This is manifested in the different types 
and numbers of fixed points generated in the renormalized dy- 
namics, and in the resulting sensitivity of the normalization 
dynamics to the initial conditions. 



VI. MANY-PARTICLE MEAN-FIELD 
CORRESPONDENCE 

Let us finally compare the mean-field description with the 
full many-particle behavior. We begin with a comparison of 
the spectral behavior. For this purpose we first have to de- 
fine the eigenenergies of the mean-field system. We will iden- 
tify them with the values of the Hamiltonian function at the 
fixed points of the mean-field dynamics. Note that these are 
different from the generalized eigenvalues of the nonlinear 
non-Hermitian Schrodinger operator (the chemical potentials) 
which were investigated in some detail for a closely related 
model in |47l l8lh . Figure [13] shows the real and imaginary 
parts of the eigenenergies in dependence on y for two different 
values of the nonlinearity. For nonvanishing nonlinearity we 
observe a similar behavior as in the linear case, where the two 
eigenvalues are purely imaginary until they meet at the criti- 
cal value |y| = v and turn into a complex conjugate pair. Here, 
however, the two eigenvalues vanish after their "collision", 
which is connected to the collision and simultaneous destruc- 
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FIG. 13: Real (left) and imaginary (right) part of the mean-field 
eigenenergies (values of the Hamiltonian function at the fixed points) 
as a function of y for 8 = 0, v = 1 and g = 0.5 (top) and g = 0.9 (bot- 
tom). 



tion of the saddle point with the center in the phase space. 
In particular, the energy values of these two fixed points are 
identical to the linear case. This is evident from the fact that 
they are located at the equator of the Bloch sphere, that is, at 
s z = and thus the nonlinear term (proportional to s\) in the 
energy vanishes. However, for values of y above the saddle- 
center collision, we still have two eigenvalues associated with 
the sink and the source that result from the bifurcatio n of one 
of the original centers at the critical value \y C nt\ = V y2 ~ 8 2 - 
Their imaginary parts are always nonzero, due to the fact that 
they are located at values s z ^ 0. Thus, the critical value for 
the emergence of the sink and the source defines the border of 
unbroken -symmetry for the mean-field system. In agree- 
ment with the many-particle results, we conclude that the non- 
linearity g shrinks the region of unbroken fPT -symmetry. 

The observed behavior is evidently the counterpart of the 
pairwise crossing structure and the unfolding of the EP of 
higher order in the many-particle spectrum. For a better 
comparison we show both the many-particle and mean-field 
eigenenergies in Fig.[l4]for the T -symmetric case as a func- 
tion of y for an intermediate interaction strength, g = 0.9. We 
indeed observe that the qualitative phenomenon of the shrink- 
ing region of unbroken !PT -symmetry is reproduced by the 
mean-field energies. However, the critical value of y that de- 
fines this border is different for the two descriptions. This 
is not surprising if we account for two facts: First, we note 
that the positions of the individual EPs depend on the particle 
number N, and the large N limit (in which one assumes the 
mean-field description to be valid) is not reached for N = 20 
particles, as in the present figure. It is in general an open ques- 
tion in which manner the mean-field limit is approached for 
non-Hermitian systems. Second, we do not expect an individ- 
ual feature of the spectrum to have an impact on the classical 
limit. This is due to the fact that this limit is only defined up 
to arbitrary orders of H (i.e. l/N in the present case), whereas 
the exact positions of individual structures, such as excep- 



tional points, is dependent on these additional terms. There- 
fore, usually isolated degeneracies do not have counterparts in 
the associated classical limit. Only if there is an accumulation 
of such points one expects a direct correspondence. Nonethe- 
less, in the present case the !PT -symmetry itself is mirrored 
in the classical system and thus we expect the breaking of this 
symmetry to be present as well. This is in agreement with 
the observed behavior for which the breaking of the symme- 
try takes place both in the mean-field and the many-particle 
system, and the influence of the interaction shrinks the region 
of unbroken !PT symmetry in both cases. 




y -2 1 

FIG. 14: (Color online) Many-particle (gray) and mean-field (dark 
red) energies for the fPT -symmetric system (8 = 0) for N = 20 par- 
ticles as functions of y, for v = 1 and g = 0.9. 

To get some insights into the correspondence of the mean- 
field and many-particle dynamics, we show several examples 
in Figs. [l5]and[T6l The figures on the top show the dynam- 
ics of the mean-field and the many-particle Bloch vector for 
a state initially located at the north pole of the Bloch sphere. 
For a better comparison we depict the dynamics of the cor- 
responding z-component, that is, the relative population im- 
balance of the two modes, in the plots in the middle. The 
resulting time dependence of the overall probability is shown 
in the lower plots. Here we have to compare the mean-field 
probability n(t), given by the normalization of the single par- 
ticle wave function, to the normalization of the many-particle 
wave function in the following way: 

n(t) = \y 1 (t)\ 2 + \y 2 \ 2 ^ a/WOW)), (50) 
thus accounting for different values of the particle number N. 

Let us first focus on Fig. Q3] where we show the dynam- 
ics for N = 20 particles. The left column shows an example 
where both the interaction strength and the non-Hermiticity 
are small (g = 0.5 and y = 0.1). The classical mean-field 
dynamics shows the typical deformed Rabi oscillations. In 
the many-particle system we observe the familiar breakdown 
behavior. Numerical results for a longer propagation sug- 
gest that the revival phenomena are strongly suppressed by 
the non-Hermiticity. The right column shows the dynam- 
ics for a stronger interaction and a stronger decay (g = 2 
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FIG. 15: (Color online) Many-particle (dashed black lines) for N = 
20 particles and mean-field (solid blue lines) dynamics for an ini- 
tial state in mode 1, that is, the north pole of the Bloch sphere with 
v = 1, e = and g = 0.5, y = 0.1 (left plots), g = 2, y = 0.5 (right 
plots). The upper plots show the dynamics of the angular momentum 
expectation values and the mean-field Bloch vector, respectively. The 
middle panels show the corresponding z-component and the lower 
plots the evolution of the overall probability n(t). 



and y = 0.5), i.e. in the mean-field self-trapping region. The 
mean-field trajectory, commenced from the north pole, ap- 
proaches the fixed point located at s z = —0.433. The full 
many-particle system shows a very similar behavior. For both 
examples the many-particle survival probability, depicted in 
the lower panel, is also reproduced by the mean-field approx- 
imation. In the regime of strong interaction, we can also ob- 
serve more complicated behavior, including phenomena re- 
lated to a many-particle tunneling from one self-trapping state 
to the other. This is illustrated in Fig. [16] where we plot the 
dynamics for large values of the interaction strength and com- 
paratively small values of y for an initial state at the north 
pole. The left column shows an example with N = 20 parti- 
cles for the parameters g = 3 and y = 0.1. One clearly ob- 
serves a tunneling of the many-particle dynamics between the 
mean-field stationary states. However, due to the fact that the 
stationary state on the south pole of the sphere is the sink of 
the mean-field dynamics, the latter approaches the southern 
fixed point as well. The right column shows a similar exam- 
ple with only N = 5 particles for a slightly smaller interaction 
strength g = 2 and a very small decay y = 0.01, to make the 
tunneling process apparent. This superimposed many-particle 
effect induces a clear mismatch into the correspondence of 
the survival-probability evolution, which is illustrated in the 
lower panels. 



FIG. 16: (Color online) Many-particle (dashed black lines) and 
mean-field (solid blue lines) dynamics for different parameters and 
an initial state in mode 1, that is, the north pole of the Bloch sphere. 
As in Fig. [l5]but for the parameters v = 1,8 = and g = 3, y = 0. 1 , 
and N = 20 particles (left plots) and g = 2, y = 0.01, and N = 5 
particles (right plots). 



The approach to the mean-field limit with increasing par- 
ticle number can be illustrated by comparison of the half life 
time of the normalization as a function of the initial condi- 
tions for different particle numbers. In Fig. [T7] we show the 
half life time as a function of the initial position on the Bloch 
sphere for y = 0.1 and g = 1 and different particle numbers. 
The corresponding mean-field behavior is depicted in the right 
plot in the middle row of Fig. Qj] with the same colorscale. It 
can be nicely seen how the mean-field features become more 
pronounced with increasing particle number. 

The presented results give a first impression on the intri- 
cate correspondence of the full many-particle description and 
the mean-field approximation for this non-Hermitian system. 
Further investigations of this correspondence and in particular 
the manner in which the mean-field limit is approached are 
promising topics for future investigations. 



VII. SUMMARY AND OUTLOOK 

We have studied the dynamics of a non-Hermitian two- 
mode Bose-Hubbard system and a related !Pf -symmetric 
model. We have derived a non-Hermitian mean-field approx- 
imation, which can be expressed in a generalized canonical 
form, including a metric gradient flow [22], and demonstrated 
the close correspondence of the damped (pseudo)classical 
motion in this mean-field description and the quantum many- 
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FIG. 17: Half life time of the rescaled normalization ^/(¥(f)|¥(f)) 
of the many-particle wave function for 8 = for v = 1, y = 0.1, 
g = l/N and different particle numbers (from left to right and top 
to bottom: N = 5, 10, 15,30). 



particle evolution. In particular, we have analyzed the fixed 
point structure of the mean-field dynamics and its bifurcation 
arising when the system parameters are varied. This results 
in a rich variety of phenomena in the many-particle dynamics, 
as for instance breakdown and revival, and tunneling, which 
can be interpreted easily in terms of the underlying mean-field 
structure. 

In conclusion, the combined presence of interaction and 
non-Hermiticity introduces a variety of interesting phenom- 
ena into the correspondence between the many-particle dy- 
namics and the mean-field description. The understanding of 
general quantum classical correspondence for non-Hermitian 
systems will ultimately require the development of new taylor- 
made methods, such as the Husimi-Schur phase space rep- 
resentation l40ll that was recently suggested in the context 
of open quantum maps. The simple model presented here 
provides an ideal testing ground for new methods for non- 
Hermitian systems. Future investigation and categorization of 
its behavior are thus a promising starting point for the for- 
mulation of a general framework for quantum classical corre- 
spondence in the presence of non-Hermiticity. 
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Appendix A: The generalized canonical equations in terms of 
the coordinates cpj 

The nonlinear complex Schrodinger equation ([3Tb can also 
be expressed in terms of a generalized canonical equation of 
motion, in the complex form 




A V#-iG" A Vr, 



(Al) 



where we have to pick two canonical conjugate variables <p 
and (p* from the four variables cpi , <p* 5 cp2 , cp| - Although the 
dynamics is apparently governed by all four variables the nor- 
malization is fixed and the dynamics is independent of the 
global phase. Therefore, we have only two independent vari- 
ables which we can choose out of the original four. It is con- 
venient to choose (pi and q>| rather than (pi and (p2- They are 
connected to the coordinates /?, q via 



9i 



1 



-2\q 



<Pl 



1 



v 2iq 



(A2) 



The equation of motion for the other variables are then implic- 
itly provided. With the choice (pi and cp* for the independent 
variables we automatically demanded (p2 to be real and ful- 
fill the normalization condition (p2 = <p2 = ~ <p*<Pi- The 
symplectic matrix is the familiar one and for the inverse of the 
Kahler metric we find: 



<Pi(l<Pil 2 -2) 2-2|( Pl | 2 +|cp 1 | 4 
2(l-|cpi| 2 ) 2(l-|cpi 
^-2|cpi| 2 +|cpi| 4 <pf(l<Pi| 2 



P) 

-2) 



. (A3) 



2(H9ll 2 ) 



2(1- 



The equations of motion for (pi and cp* can then be found from 
(lAlb . where H and T are given by the real and imaginary parts 
of the Hamiltonian function for the non-Hermitian and nonlin- 
ear two-level system expressed in terms of (pi and (pj: 



^ = (e-iY)(cpIcp 1 -l)+v v /l-cp*cpi(cpt + cpi)+|(cptcp 1 -l) 2 . 

(A4) 

The equation of motion for 92 can then be deduced from the 
dynamics of (pi via 



q>2 



(j)i(p* + (pi(p* 
2 v /1 -9i9t 



(A5) 



The dynamics thus obtained is equivalent to the non- 
Hermitian (nonlinear) Schrodinger equation (l37l) up to a 
global phase. 
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